Multinomial Logistic Regression

Advanced Categorical Data Analysis

Author

Dr Muhammad Saufi

Published

July 18, 2024

This document serves both learning and practical purposes. It is designed for educational use, aiming to enhance statistical analysis skills and provide clear, organized notes for future reference.

1 Introduction

Multinomial Logistic Regression is an extension of binary logistic regression that allows for the modeling of outcome variables with more than two categories. This type of regression is used when the outcome categories are nominal, meaning they do not have a natural order. The dependent variable has multiple categories (more than two), and one of the categories is chosen as the reference category. The model compares each of the other categories to this reference category.

1.1 Model Structure

  • The model estimates the log odds of each category relative to the reference category.
  • For a categorical outcome with \(k\) categories, the model fits \(k-1\) equations.
  • Each equation models the log odds of a specific category relative to the reference category.

1.2 Model Formulation

  1. Log Odds
  • The log odds of being in category \(j\) relative to the reference category \(0\) is modeled as:

\[ \log \left( \frac{P(Y=j|X)}{P(Y=0|X)} \right) = \alpha_j + \beta_{j1}X_1 + \beta_{j2}X_2 + \ldots + \beta_{jp}X_p \]

Note

\(\alpha_j\) is the intercept for category \(j\), and \(\beta_{jk}\) are the coefficients for the predictor variables \(X_k\).

  1. Probability
  • The probability of being in category \(j\) is derived from the log odds:

\[ P(Y=j|X) = \frac{e^{\alpha_j + \beta_{j1}X_1 + \beta_{j2}X_2 + \ldots + \beta_{jp}X_p}}{1 + \sum_{l=1}^{k-1} e^{\alpha_l + \beta_{l1}X_1 + \beta_{l2}X_2 + \ldots + \beta_{lp}X_p}} \]

  • The probability of being in the reference category \(0\) is:

\[ P(Y=0|X) = \frac{1}{1 + \sum_{l=1}^{k-1} e^{\alpha_l + \beta_{l1}X_1 + \beta_{l2}X_2 + \ldots + \beta_{lp}X_p}} \]

1.3 Estimation

  • The parameters (\(\alpha_j\) and \(\beta_{jk}\)) are estimated using maximum likelihood estimation (MLE).
  • The likelihood function is constructed based on the observed data and the probabilities of each category.
  • Iterative algorithms, such as Newton-Raphson or Fisher scoring, are used to find the parameter estimates that maximize the likelihood function.

1.4 Interpretation

  1. Coefficients (\(\beta\))

    • The coefficients represent the change in the log odds of being in a category relative to the reference category for a one-unit change in the predictor variable.
    • Positive coefficients indicate higher odds of being in the category compared to the reference, while negative coefficients indicate lower odds.
  2. Odds Ratios (OR) is obtained by exponentiating the coefficients:

\[ OR = e^{\beta} \]

Note

An OR greater than 1 indicates increased odds of being in the category compared to the reference, and an OR less than 1 indicates decreased odds.

1.5 Model Fit and Diagnostics

  • Likelihood Ratio Test: Used to compare nested models and assess the significance of adding or removing predictors.
  • Wald Test: Used to test the significance of individual coefficients.

2 Setting Up the Environment

Load the necessary libraries to ensure the R environment is equipped with the tools and functions required for efficient and effective analysis.

3 Practical 1: Mammography Experience Dataset

In this exercise, the objectives are to understand and apply multinomial logistic regression using a dataset related to mammography experiences. The focus is on analyzing the mammography experiences dataset to understand factors influencing whether women have had a mammogram recently, a long time ago, or never.

3.1 Dataset Overview

The dataset used is focused on mammography experiences which consists of responses from a survey on mammography experiences and perceptions. The dataset includes the following variables:

  1. me (Mammography Experience): This is the outcome variable and is categorical with three levels:

    • 0: Never had a mammogram
    • 1: Had a mammogram within a year
    • 2: Had a mammogram over a year ago
  2. sympt (Symptoms Perception): This variable captures the perception that a mammogram is not needed unless symptoms are present.

    • 1: Strongly Agree
    • 2: Agree
    • 3: Disagree
    • 4: Strongly Disagree
  3. pb (Perceived Benefit of Mammography): This variable measures the perceived benefit of mammography on a scale. The exact scale is not specified but typically ranges from lower to higher perceived benefits.

  4. hist (Family History of Breast Cancer): This binary variable indicates whether the respondent has a family history (mother or sister) of breast cancer.

    • 0: No
    • 1: Yes
  5. bse (Breast Self-Examination Training): This binary variable indicates whether the respondent has been taught how to perform breast self-examinations.

    • 0: No
    • 1: Yes
  6. detc (Detection Likelihood): This variable measures the respondent’s perception of the likelihood that a mammogram could detect a new case of breast cancer. It is coded as:

    • 1: Not likely
    • 2: Somewhat likely
    • 3: Very likely

3.2 Data Reading and Exploration

# Read the data
dat.m1 <- read_dta("Datasets/mammog9.dta")

# Explore the data
glimpse(dat.m1)
Rows: 412
Columns: 7
$ obs   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
$ me    <dbl> 0, 0, 0, 1, 2, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0…
$ sympt <dbl> 3, 2, 3, 3, 4, 3, 4, 4, 2, 4, 4, 3, 4, 1, 2, 4, 3, 4, 3, 3, 3, 4…
$ pb    <dbl> 7, 11, 8, 11, 7, 7, 6, 6, 6, 6, 8, 6, 6, 5, 8, 11, 6, 5, 10, 10,…
$ hist  <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
$ bse   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1…
$ detc  <dbl> 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 2…
summary(dat.m1)
      obs              me             sympt             pb        
 Min.   :  1.0   Min.   :0.0000   Min.   :1.000   Min.   : 5.000  
 1st Qu.:103.8   1st Qu.:0.0000   1st Qu.:2.000   1st Qu.: 6.000  
 Median :206.5   Median :0.0000   Median :3.000   Median : 7.000  
 Mean   :206.5   Mean   :0.6117   Mean   :2.966   Mean   : 7.561  
 3rd Qu.:309.2   3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.: 9.000  
 Max.   :412.0   Max.   :2.0000   Max.   :4.000   Max.   :17.000  
      hist             bse              detc      
 Min.   :0.0000   Min.   :0.0000   Min.   :1.000  
 1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:2.000  
 Median :0.0000   Median :1.0000   Median :3.000  
 Mean   :0.1068   Mean   :0.8689   Mean   :2.658  
 3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:3.000  
 Max.   :1.0000   Max.   :1.0000   Max.   :3.000  

3.3 Data Wrangling

The me and hist variables are treated as numerical despite being categorical variables. To ensure that the data is handled correctly during statistical analyses, these variables must be converted to factors so that each level is treated as a distinct category, and the analysis methods can correctly interpret these categories.

3.3.1 me Variable

# Counting the Original `me` Variable
dat.m1 %>% count(me)
# A tibble: 3 × 2
     me     n
  <dbl> <int>
1     0   234
2     1   104
3     2    74
# Copy dataset and keep the original dataset
dat.m2 <- dat.m1

# Converting `me` to a factor
dat.m2 <- dat.m2 %>%
  mutate(me2 = factor(me, labels = c("never", "within.a.year", "over.a.year.ago")))
dat.m2 %>% count(me2)
# A tibble: 3 × 2
  me2                 n
  <fct>           <int>
1 never             234
2 within.a.year     104
3 over.a.year.ago    74
Note

This code creates a new variable me2 by converting the numeric me variable to a factor. The labels argument assigns meaningful labels to the factor levels.

3.3.2 hist Variable

# Counting the Original `hist` Variable
dat.m2 %>% count(hist)
# A tibble: 2 × 2
   hist     n
  <dbl> <int>
1     0   368
2     1    44
# Converting `hist` to a Factor
dat.m2 <- dat.m2 %>% mutate(hist2 = factor(hist, labels = c("no", "yes")))
dat.m2 %>% count(hist2)
# A tibble: 2 × 2
  hist2     n
  <fct> <int>
1 no      368
2 yes      44
Note

This code creates a new variable hist2 by converting the numeric hist variable to a factor, with no and yes as labels.

3.3.3 Table Summary

# Generate summary statistics for the dataset, grouped by the `me2` variable
dat.m2 %>%
  select(-me, -hist) %>%
  tbl_summary(by = me2, statistic = list(all_continuous() ~ "{mean} ({sd})"))
Characteristic never, N = 2341 within.a.year, N = 1041 over.a.year.ago, N = 741
obs 203 (121) 212 (113) 211 (122)
sympt


    1 33 (14%) 2 (1.9%) 5 (6.8%)
    2 62 (26%) 4 (3.8%) 7 (9.5%)
    3 85 (36%) 43 (41%) 32 (43%)
    4 54 (23%) 55 (53%) 30 (41%)
pb 8.06 (2.19) 6.70 (1.66) 7.19 (1.91)
bse 190 (81%) 99 (95%) 69 (93%)
detc


    1 13 (5.6%) 1 (1.0%) 4 (5.4%)
    2 77 (33%) 12 (12%) 16 (22%)
    3 144 (62%) 91 (88%) 54 (73%)
hist2 14 (6.0%) 19 (18%) 11 (15%)
1 Mean (SD); n (%)

3.4 Estimation

The VGAM package is used to run multinomial logistic regression, which is suitable for a dependent variable with more than two categories. The objctive is to use the VGAM package to fit a multinomial logistic regression model and ensure the reference category for the dependent variable me2 is “never”.

3.4.1 Confirm the Order of Factor Levels

Check the current order of the factor levels for me2. It is crucial because the reference level in regression models affects the interpretation of coefficients.

levels(dat.m2$me2)
[1] "never"           "within.a.year"   "over.a.year.ago"
Note

The levels are currently in the order where “never” is the first level. For the model, we need “never” to be the last level to serve as the reference category.

3.4.2 Reorder the Factor Levels

# Create new dataset `dat.m3`
dat.m3 <- dat.m2

# Generate `me3`
dat.m3 <- dat.m3 %>%
  mutate(me3 = fct_relevel(me2, c("over.a.year.ago", "within.a.year", "never")))
Note

Creates a new variable me3 by reordering the levels of me2 using fct_relevel so that “never” becomes the reference level.

table(dat.m3$me)

  0   1   2 
234 104  74 

summary(dat.m3$me2)
          never   within.a.year over.a.year.ago 
            234             104              74 

summary(dat.m3$me3)
over.a.year.ago   within.a.year           never 
             74             104             234 
Note

These functions confirm that the observations and levels are correctly adjusted after releveling.

3.4.3 Check the Order of Levels

levels(dat.m3$me2)
[1] "never"           "within.a.year"   "over.a.year.ago"

levels(dat.m3$me3)
[1] "over.a.year.ago" "within.a.year"   "never"          
Note

This checks the new order of levels to ensure that me3 has “never” as the last level (reference category).

3.4.4 Multinomial Logistic Regression Model

3.4.4.1 Model with hist2 as Predictor

fitmlog1 <- vglm(me3 ~ hist2,
  family = multinomial(),
  data = dat.m3
)
summary(fitmlog1)

Call:
vglm(formula = me3 ~ hist2, family = multinomial(), data = dat.m3)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  -1.2505     0.1429  -8.751  < 2e-16 ***
(Intercept):2  -0.9510     0.1277  -7.446  9.6e-14 ***
hist2yes:1      1.0093     0.4275   2.361 0.018225 *  
hist2yes:2      1.2564     0.3747   3.353 0.000798 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 792.3399 on 820 degrees of freedom

Log-likelihood: -396.17 on 820 degrees of freedom

Number of Fisher scoring iterations: 4 

No Hauck-Donner effect found in any of the estimates


Reference group is level  3  of the response

  • vglm(me3 ~ hist2, family = multinomial(), data = dat.m1): Fits a model where me3 (mammography experience) is the dependent variable, and hist2 (family history of breast cancer) is the predictor.
  • family = multinomial(): Specifies that the model is a multinomial logistic regression.
  • data = dat.m3: Specifies the dataset to be used.

Summary of Output

  1. Intercepts

    • (Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” when hist2 is the reference level (no family history).
    • (Intercept):2: The log-odds of being in the “within a year” category relative to “never” when hist2 is the reference level (no family history).
  2. Predictor (hist2yes)

    • hist2yes:1: The log-odds change of being in the “over a year ago” category relative to “never” for participants with a family history of breast cancer compared to those without.
    • hist2yes:2: The log-odds change of being in the “within a year” category relative to “never” for participants with a family history of breast cancer compared to those without.
  3. Statistical Significance

    • (Intercept):1 and (Intercept):2 are highly significant (p < 0.001), indicating that the intercepts are significantly different from zero.
    • hist2yes:1 is significant at the 0.05 level (*), indicating that the log-odds change for the “over a year ago” category relative to “never” is significant for those with a family history of breast cancer.
    • hist2yes:2 is highly significant (p < 0.001), indicating that the log-odds change for the “within a year” category relative to “never” is significant for those with a family history of breast cancer.
  4. Interpretation

    • The model indicates that having a family history of breast cancer (hist2 = “yes”) significantly affects the likelihood of having had a mammogram “within a year” or “over a year ago” compared to “never”.
    • Significant positive coefficients for hist2yes suggest that participants with a family history of breast cancer are more likely to have had a mammogram within a year or over a year ago compared to those without a family history.
    • The model fit statistics (residual deviance, log-likelihood) suggest that the model adequately fits the data.

3.4.4.2 Model with detc as Predictor

# Convert `detc` to a factor
dat.m3 <- dat.m3 %>% mutate(detc2 = factor(detc))

# Check the distribution of `detc2` across the levels of `me3`
table(dat.m3$me3, dat.m3$detc2)
                 
                    1   2   3
  over.a.year.ago   4  16  54
  within.a.year     1  12  91
  never            13  77 144
Note

The table shows the count of observations for each combination of me2 and detc2.

# Fitting the model with `detc2` as predictor
fitmlog2 <- vglm(me3 ~ detc2,
  family = multinomial(),
  data = dat.m3
)
summary(fitmlog2)

Call:
vglm(formula = me3 ~ detc2, family = multinomial(), data = dat.m3)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)  
(Intercept):1  -1.1787     0.5718  -2.061   0.0393 *
(Intercept):2  -2.5649     1.0377      NA       NA  
detc22:1       -0.3926     0.6344  -0.619   0.5360  
detc22:2        0.7061     1.0832   0.652   0.5145  
detc23:1        0.1978     0.5936   0.333   0.7389  
detc23:2        2.1060     1.0463   2.013   0.0441 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 778.4011 on 818 degrees of freedom

Log-likelihood: -389.2005 on 818 degrees of freedom

Number of Fisher scoring iterations: 5 

Warning: Hauck-Donner effect detected in the following estimate(s):
'(Intercept):2'


Reference group is level  3  of the response

Summary of Output

  1. Intercepts

    • (Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” when detc2 is at its reference level (usually the first level).
    • (Intercept):2: The log-odds of being in the “within a year” category relative to “never” when detc2 is at its reference level.
  2. Predictor (detc2)

    • detc22:1: The log-odds change of being in the “over a year ago” category relative to “never” for participants with the second level of detc2 compared to the reference level.
    • detc22:2: The log-odds change of being in the “within a year” category relative to “never” for participants with the second level of detc2 compared to the reference level.
    • detc23:1: The log-odds change of being in the “over a year ago” category relative to “never” for participants with the third level of detc2 compared to the reference level.
    • detc23:2: The log-odds change of being in the “within a year” category relative to “never” for participants with the third level of detc2 compared to the reference level.
  3. **Statistical Significance*

    • (Intercept):1 is significant at the 0.05 level (*), indicating that the intercept for the “over a year ago” category relative to “never” is significantly different from zero.
    • detc23:2 is significant at the 0.05 level (*), indicating that the log-odds change for the “within a year” category relative to “never” is significant for those at the third level of detc2.
  4. Interpretation

    • The model indicates that the perceived likelihood of detection (detc2) has an effect on the likelihood of having had a mammogram “within a year” compared to “never”.
    • Significant positive coefficient for detc23:2 suggests that participants with the highest level of perceived likelihood of detection are more likely to have had a mammogram within a year compared to those with the lowest level of perceived likelihood.
    • The model fit statistics (residual deviance, log-likelihood) suggest that the model adequately fits the data.

3.5 Inferences

Inferences refer to the process of drawing conclusions about a population based on the analysis of sample data. In multinomial logistic regression model, inferences involve estimating the relationships between the predictor variables and the outcome variable and determining the statistical significance and confidence of these estimates. The process of inferences involves:

  1. Calculate the 95% Confidence Intervals (CIs) for the model coefficients.
  2. Calculate the p-values for hypothesis testing.
  3. Exponentiate the coefficients to obtain the Relative Risk Ratios (RRRs).

3.5.1 Calculate 95% Confidence Intervals and Combine with Coefficients

# Extract coefficients from the model
b_fitmlog2 <- coef(fitmlog2)

# Calculate confidence intervals
ci_fitmlog2 <- confint(fitmlog2)

# Combine coefficients and confidence intervals
b_ci_fitmlog2 <- cbind(b_fitmlog2, ci_fitmlog2)
b_ci_fitmlog2
              b_fitmlog2       2.5 %      97.5 %
(Intercept):1 -1.1786550 -2.29930795 -0.05800204
(Intercept):2 -2.5649494 -4.59888160 -0.53101711
detc22:1      -0.3925617 -1.63588117  0.85075776
detc22:2       0.7060506 -1.41689336  2.82899454
detc23:1       0.1978257 -0.96565093  1.36130242
detc23:2       2.1059956  0.05519791  4.15679321

  • coef(fitmlog2): Extracts the estimated coefficients from the model fitmlog2.
  • confint(fitmlog2): Calculates the 95% confidence intervals for the model coefficients.
  • cbind(b_fitmlog2, ci_fitmlog2): Combines the coefficients and their confidence intervals into a single table.
  • b_fitmlog2: The estimated coefficients.
  • 2.5 % and 97.5 %: The lower and upper bounds of the 95% confidence intervals for each coefficient.

3.5.2 Exponentiate Coefficients to Obtain Relative Risk Ratios (RRRs)

# Exponentiate the combined coefficients and confidence intervals
rrr_fitmlog2 <- exp(b_ci_fitmlog2)

# Combine the RRRs with the original table
tab_fitmlog2 <- cbind(b_ci_fitmlog2, rrr_fitmlog2)

# Name the columns of the combined table
colnames(tab_fitmlog2) <- c("b", "lower b", "upper b", "rrr", "lower rrr", "upper rrr")

# Display the combined table
tab_fitmlog2
                       b     lower b     upper b        rrr  lower rrr
(Intercept):1 -1.1786550 -2.29930795 -0.05800204 0.30769231 0.10032825
(Intercept):2 -2.5649494 -4.59888160 -0.53101711 0.07692308 0.01006308
detc22:1      -0.3925617 -1.63588117  0.85075776 0.67532468 0.19478066
detc22:2       0.7060506 -1.41689336  2.82899454 2.02597403 0.24246610
detc23:1       0.1978257 -0.96565093  1.36130242 1.21875000 0.38073529
detc23:2       2.1059956  0.05519791  4.15679321 8.21527778 1.05674974
               upper rrr
(Intercept):1  0.9436480
(Intercept):2  0.5880066
detc22:1       2.3414204
detc22:2      16.9284313
detc23:1       3.9012711
detc23:2      63.8663880

  • exp(b_ci_fitmlog2): Exponentiates the coefficients and their confidence intervals to obtain the RRRs.
  • cbind(b_ci_fitmlog2, rrr_fitmlog2): Combines the original table with the RRRs.
  • colnames(tab_fitmlog2) <- c('b', 'lower b', 'upper b', 'rrr', 'lower rrr', 'upper rrr'): Renames the columns of the combined table for clarity.
  • b: The estimated coefficients.
  • lower b and upper b: The lower and upper bounds of the 95% confidence intervals for the coefficients.
  • rrr: The Relative Risk Ratios (exponentiated coefficients).
  • lower rrr and upper rrr: The lower and upper bounds of the 95% confidence intervals for the RRRs.

Interpretation


  • (Intercept):1

    • The log-odds of being in the “over a year ago” category relative to “never” are -1.1786550 when all predictors are at their reference levels.
    • The confidence interval for this log-odds is [-2.29930795, -0.0580204], indicating it is statistically significant since it does not include zero.
    • The relative risk ratio (RRR) is 0.3076921, meaning that the likelihood of being in the “over a year ago” category is about 0.31 times that of the “never” category.
    • The 95% confidence interval for the RRR is [0.10032825, 0.9434680], indicating the precision of the estimate.
  • detc22:1

    • The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.3925617 for each one-unit increase in detc22.
    • The confidence interval for this coefficient is [-1.63588117, 0.85075776], which includes zero, indicating it is not statistically significant.
    • The RRR is 0.67532468, meaning the likelihood decreases to about 0.68 times for each one-unit increase in detc22.
    • The 95% confidence interval for the RRR is [0.19478066, 2.3414204].
  • detc23:1

    • The log-odds of being in the “over a year ago” category relative to “never” increase by 0.1978257 for each one-unit increase in detc23.
    • The confidence interval for this coefficient is [-0.96556939, 1.36130242], which includes zero, indicating it is not statistically significant.
    • The RRR is 1.21875000, meaning the likelihood increases to about 1.22 times for each one-unit increase in detc23.
    • The 95% confidence interval for the RRR is [0.38073529, 3.9012711].

  • (Intercept):2

    • The log-odds of being in the “within a year” category relative to “never” are -2.5649494 when all predictors are at their reference levels.
    • The confidence interval for this log-odds is [-4.59888160, -0.5310171], indicating it is statistically significant.
    • The RRR is 0.07692308, meaning that the likelihood of being in the “within a year” category is about 0.08 times that of the “never” category.
    • The 95% confidence interval for the RRR is [0.01006308, 0.5888066].
  • detc22:2

    • The log-odds of being in the “within a year” category relative to “never” increase by 0.7060506 for each one-unit increase in detc22.
    • The confidence interval for this coefficient is [-1.41689336, 2.82899454], which includes zero, indicating it is not statistically significant.
    • The RRR is 2.02597403, meaning the likelihood increases to about 2.03 times for each one-unit increase in detc22.
    • The 95% confidence interval for the RRR is [0.24246610, 16.9284133].
  • detc23:2

    • The log-odds of being in the “within a year” category relative to “never” increase by 2.1059956 for each one-unit increase in detc23.
    • The confidence interval for this coefficient is [0.05519791, 4.15679321], indicating it is statistically significant since it does not include zero.
    • The RRR is 8.21527778, meaning the likelihood increases to about 8.22 times for each one-unit increase in detc23.
    • The 95% confidence interval for the RRR is [1.05674974, 63.8663880], indicating a high degree of variability in the estimate.

4 Practical 2: Additional Covariates

In this exercise, the objective is to extend the initial multinomial logistic regression model by adding more covariates to better understand the relationships between the predictors and the outcome variable. This helps in refining the model and potentially improving its explanatory power.

4.1 Estimation

Define the additional covariates to be included in the model. These include variables such as:

  • sympt (perception of symptoms)
  • pb (perceived benefit)
  • hist (family history)
  • bse (breast self-exam training)
  • detc (detection likelihood)

4.1.1 Defining the Model

The model is specified as a multinomial logistic regression, with the dataset dat.m3. The additional covariates are included in the model.

# Fitting the model with additional covariates
fitmlog3 <- vglm(me3 ~ factor(sympt) + pb + hist2 + bse + detc2,
  family = multinomial(),
  data = dat.m3
)
summary(fitmlog3)

Call:
vglm(formula = me3 ~ factor(sympt) + pb + hist2 + bse + detc2, 
    family = multinomial(), data = dat.m3)

Coefficients: 
                 Estimate Std. Error z value Pr(>|z|)   
(Intercept):1    -0.98609    1.11184  -0.887  0.37513   
(Intercept):2    -2.99875    1.53905      NA       NA   
factor(sympt)2:1 -0.29008    0.64406  -0.450  0.65243   
factor(sympt)2:2  0.11004    0.92273   0.119  0.90508   
factor(sympt)3:1  0.81731    0.53979   1.514  0.12999   
factor(sympt)3:2  1.92471    0.77757   2.475  0.01331 * 
factor(sympt)4:1  1.13224    0.54767   2.067  0.03870 * 
factor(sympt)4:2  2.45699    0.77530   3.169  0.00153 **
pb:1             -0.14821    0.07637  -1.941  0.05230 . 
pb:2             -0.21944    0.07551  -2.906  0.00366 **
hist2yes:1        1.06544    0.45940   2.319  0.02038 * 
hist2yes:2        1.36624    0.43752   3.123  0.00179 **
bse:1             1.05214    0.51499   2.043  0.04105 * 
bse:2             1.29167    0.52988   2.438  0.01478 * 
detc22:1         -0.92439    0.71375  -1.295  0.19528   
detc22:2          0.01702    1.16169   0.015  0.98831   
detc23:1         -0.69053    0.68712  -1.005  0.31491   
detc23:2          0.90414    1.12661   0.803  0.42225   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 693.9019 on 806 degrees of freedom

Log-likelihood: -346.951 on 806 degrees of freedom

Number of Fisher scoring iterations: 5 

Warning: Hauck-Donner effect detected in the following estimate(s):
'(Intercept):2'


Reference group is level  3  of the response

  • vglm(): Fits a multinomial logistic regression model with me3 as the dependent variable and the specified covariates as predictors.
  • factor(sympt) converts sympt to a factor variable.

4.1.2 Checking and Recoding sympt Variable

The variable sympt has 4 categories, hence it is needed to be recoded into a binary factor variable symptd for simplicity.

# Original `sympt` variable
dat.m3 %>% count(sympt)
# A tibble: 4 × 2
  sympt     n
  <dbl> <int>
1     1    40
2     2    73
3     3   160
4     4   139
# recode `sympt` into a binary variable `symptd`
dat.m3 <- dat.m3 %>%
  mutate(symptd = ifelse(sympt < 3, "<3", "3 or more"))
dat.m3 %>% count(sympt, symptd)
# A tibble: 4 × 3
  sympt symptd        n
  <dbl> <chr>     <int>
1     1 <3           40
2     2 <3           73
3     3 3 or more   160
4     4 3 or more   139

  • mutate(symptd = ifelse(sympt < 3, '<3', '3 or more')): Recode sympt into a binary variable symptd, where categories less than 3 are labeled <3 and categories 3 or more are labeled 3 or more.

4.1.3 Refit the Model with New symptd Variable

# Refit the model using the new binary variable `symptd`
fitmlog4 <- vglm(me3 ~ symptd + pb + hist2 + bse + detc2,
  family = multinomial(),
  data = dat.m3
)
summary(fitmlog4)

Call:
vglm(formula = me3 ~ symptd + pb + hist2 + bse + detc2, family = multinomial(), 
    data = dat.m3)

Coefficients: 
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept):1     -0.99877    1.07197  -0.932 0.351485    
(Intercept):2     -2.70375    1.43422  -1.885 0.059406 .  
symptd3 or more:1  1.12136    0.35720   3.139 0.001693 ** 
symptd3 or more:2  2.09534    0.45739   4.581 4.62e-06 ***
pb:1              -0.16811    0.07417  -2.266 0.023426 *  
pb:2              -0.25101    0.07293  -3.442 0.000578 ***
hist2yes:1         1.01406    0.45381   2.235 0.025446 *  
hist2yes:2         1.29328    0.43353   2.983 0.002853 ** 
bse:1              1.02859    0.51398   2.001 0.045366 *  
bse:2              1.24397    0.52630   2.364 0.018097 *  
detc22:1          -0.90213    0.71463  -1.262 0.206813    
detc22:2           0.09028    1.16079   0.078 0.938011    
detc23:1          -0.66982    0.68759  -0.974 0.329979    
detc23:2           0.97281    1.12604   0.864 0.387627    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 697.4959 on 810 degrees of freedom

Log-likelihood: -348.748 on 810 degrees of freedom

Number of Fisher scoring iterations: 5 

No Hauck-Donner effect found in any of the estimates


Reference group is level  3  of the response

Summary of Output

  1. Intercepts

    • (Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” are -0.99877 when all predictors are at their reference levels. This is not statistically significant (p = 0.351485).
    • (Intercept):2: The log-odds of being in the “within a year” category relative to “never” are -2.70375 when all predictors are at their reference levels. This is marginally significant (p = 0.059406).
  2. Predictors

    • symptd3 or more:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.12136 for those with symptoms 3 or more compared to less than 3. This is statistically significant (p = 0.001693).

    • symptd3 or more:2: The log-odds of being in the “within a year” category relative to “never” increase by 2.09534 for those with symptoms 3 or more compared to less than 3. This is highly significant (p < 0.001).

    • pb:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.16811 for each unit increase in perceived benefit. This is statistically significant (p = 0.023426).

    • pb:2: The log-odds of being in the “within a year” category relative to “never” decrease by 0.25101 for each unit increase in perceived benefit. This is highly significant (p < 0.001).

    • hist2yes:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.01406 for those with a family history of breast cancer compared to those without. This is statistically significant (p = 0.025446).

    • hist2yes:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.29383 for those with a family history of breast cancer compared to those without. This is highly significant (p = 0.002835).

    • bse:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.28509 for those with breast self-exam training compared to those without. This is statistically significant (p = 0.012473).

    • bse:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.24397 for those with breast self-exam training compared to those without. This is statistically significant (p = 0.018097).

    • detc22:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.90213 for those in the second level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.206813).

    • detc22:2: The log-odds of being in the “within a year” category relative to “never” increase by 0.09028 for those in the second level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.938011).

    • detc23:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.66982 for those in the third level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.329979).

    • detc23:2: The log-odds of being in the “within a year” category relative to “never” increase by 0.97281 for those in the third level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.387627).

  3. Interpretation

    • The model includes additional covariates such as symptoms, perceived benefit, family history, breast self-exam training, and detection likelihood.
    • Significant predictors include symptd, pb, hist2, and bse, indicating their impact on the likelihood of having a mammogram within a year or over a year ago relative to never.
    • The model fit statistics suggest that the model adequately fits the data.
    • Confidence intervals and RRRs provide further insights into the effects of these predictors.

4.1.4 Convert detc to a Binary Variable

The variable detc initially has 3 categories. The aim is to recode it into a binary variable detcd.

dat.m3 <- dat.m3 %>%
  mutate(detcd = ifelse(detc < 3, "<3", "3 or more"))

  • ifelse(detc < 3, '<3', '3 or more'): This function checks the value of detc for each observation.
  • If detc is less than 3, it assigns the value ‘<3’ to detcd.
  • If detc is 3 or more, it assigns the value ‘3 or more’ to detcd.
  • This effectively recodes detc from 3 categories into a binary variable with 2 categories: ‘<3’ and ‘3 or more’.

4.1.5 Refit the Model with the New detcd Variable

fitmlog5 <- vglm(me3 ~ symptd + pb + hist2 + bse + detcd,
  family = multinomial(),
  data = dat.m3
)
summary(fitmlog5)

Call:
vglm(formula = me3 ~ symptd + pb + hist2 + bse + detcd, family = multinomial(), 
    data = dat.m3)

Coefficients: 
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept):1     -1.82388    0.85509  -2.133 0.032928 *  
(Intercept):2     -2.62376    0.92639  -2.832 0.004622 ** 
symptd3 or more:1  1.12742    0.35636   3.164 0.001558 ** 
symptd3 or more:2  2.09475    0.45742   4.579 4.66e-06 ***
pb:1              -0.15432    0.07262  -2.125 0.033587 *  
pb:2              -0.24947    0.07258  -3.437 0.000588 ***
hist2yes:1         1.06318    0.45284   2.348 0.018885 *  
hist2yes:2         1.30986    0.43360   3.021 0.002520 ** 
bse:1              0.95601    0.50734   1.884 0.059515 .  
bse:2              1.23701    0.52542   2.354 0.018557 *  
detcd3 or more:1   0.11416    0.31821   0.359 0.719786    
detcd3 or more:2   0.88518    0.35624   2.485 0.012962 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])

Residual deviance: 699.1326 on 812 degrees of freedom

Log-likelihood: -349.5663 on 812 degrees of freedom

Number of Fisher scoring iterations: 5 

No Hauck-Donner effect found in any of the estimates


Reference group is level  3  of the response

Summary of Output

  1. Intercepts

    • (Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” are -1.82388 when all predictors are at their reference levels. This is statistically significant (p = 0.032928).
    • (Intercept):2: The log-odds of being in the “within a year” category relative to “never” are -2.62376 when all predictors are at their reference levels. This is highly significant (p = 0.004622).
  2. Predictors

    • symptd3 or more:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.12742 for those with symptoms 3 or more compared to less than 3. This is statistically significant (p = 0.001558).

    • symptd3 or more:2: The log-odds of being in the “within a year” category relative to “never” increase by 2.09475 for those with symptoms 3 or more compared to less than 3. This is highly significant (p < 0.001).

    • pb:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.15432 for each unit increase in perceived benefit. This is statistically significant (p = 0.033587).

    • pb:2: The log-odds of being in the “within a year” category relative to “never” decrease by 0.24947 for each unit increase in perceived benefit. This is highly significant (p < 0.001).

    • hist2yes:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.06318 for those with a family history of breast cancer compared to those without. This is statistically significant (p = 0.018885).

    • hist2yes:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.38968 for those with a family history of breast cancer compared to those without. This is highly significant (p = 0.001520).

    • bse:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 0.95601 for those with breast self-exam training compared to those without. This is marginally significant (p = 0.059515).

    • bse:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.20137 for those with breast self-exam training compared to those without. This is statistically significant (p = 0.018557).

    • detcd3 or more:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 0.11416 for those with detection likelihood 3 or more compared to less than 3. This is not statistically significant (p = 0.719786).

    • detcd3 or more:2: The log-odds of being in the “within a year” category relative to “never” increase by 0.88518 for those with detection likelihood 3 or more compared to less than 3. This is statistically significant (p = 0.012962).

  3. Interpretation

    • The model includes covariates such as symptoms, perceived benefit, family history, breast self-exam training, and detection likelihood.
    • Significant predictors include symptd, pb, hist2, and bse, indicating their impact on the likelihood of having a mammogram within a year or over a year ago relative to never.
    • The model fit statistics suggest that the model adequately fits the data.
    • Confidence intervals and RRRs provide further insights into the effects of these predictors.

4.2 Inferences

  1. Retrieve the log-odds and RRRs (Relative Risk Ratios).
  2. Calculate the 95% Confidence Intervals (CIs) for both the log-odds and the RRRs.
  3. Combine these results into a table.
  4. Present the results in a table

4.2.1 Retrieve Log-Odds and RRRs

The log-odds are the coefficients from the multinomial logistic regression model, and the RRRs are the exponentiated coefficients.

b_rrr_fitmlog5 <- cbind(coef(fitmlog5), exp(coef(fitmlog5)))
  • coef(fitmlog5): Extracts the coefficients (log-odds) from the model fitmlog5.
  • exp(coef(fitmlog5)): Exponentiates the coefficients to obtain the RRRs.
  • cbind(coef(fitmlog5), exp(coef(fitmlog5))): Combines the log-odds and RRRs into one table.

4.2.2 Calculate 95% Confidence Intervals

ci_b_rrr_fitmlog5 <- cbind(confint(fitmlog5), exp(confint(fitmlog5)))
  • confint(fitmlog5): Computes the 95% CIs for the log-odds.
  • exp(confint(fitmlog5)): Computes the 95% CIs for the RRRs by exponentiating the CIs for the log-odds.
  • cbind(confint(fitmlog5), exp(confint(fitmlog5))): Combines the CIs for log-odds and RRRs into one table.

4.2.3 Combine the Results into a Table

res_fitmlog5 <- cbind(b_rrr_fitmlog5, ci_b_rrr_fitmlog5)
colnames(res_fitmlog5) <- c("b", "rrr", "lower 95% b", "upper 95% b", "lower 95% rrr", "upper 95% rrr")
  • cbind(b_rrr_fitmlog5, ci_b_rrr_fitmlog5): Combines the log-odds, RRRs, and their CIs.
  • colnames(res_fitmlog5) <- c('b', 'rrr', 'lower 95% b', 'upper 95% b', 'lower 95% rrr', 'upper 95% rrr'): Renames the columns of the combined table for clarity.

4.2.4 Display the Results

kable(res_fitmlog5, digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
b rrr lower 95% b upper 95% b lower 95% rrr upper 95% rrr
(Intercept):1 -1.824 0.161 -3.500 -0.148 0.030 0.862
(Intercept):2 -2.624 0.073 -4.439 -0.808 0.012 0.446
symptd3 or more:1 1.127 3.088 0.429 1.826 1.536 6.208
symptd3 or more:2 2.095 8.123 1.198 2.991 3.314 19.911
pb:1 -0.154 0.857 -0.297 -0.012 0.743 0.988
pb:2 -0.249 0.779 -0.392 -0.107 0.676 0.898
hist2yes:1 1.063 2.896 0.176 1.951 1.192 7.034
hist2yes:2 1.310 3.706 0.460 2.160 1.584 8.669
bse:1 0.956 2.601 -0.038 1.950 0.962 7.031
bse:2 1.237 3.445 0.207 2.267 1.230 9.649
detcd3 or more:1 0.114 1.121 -0.510 0.738 0.601 2.091
detcd3 or more:2 0.885 2.423 0.187 1.583 1.206 4.871

  • library(knitr) and library(kableExtra): Load the necessary libraries for table formatting.
  • kable(res_fitmlog5, digits = 3): Creates a table with res_fitmlog5 data, formatting the values to three decimal places.
  • kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")): Applies styling options to the table for better presentation.

Interpretation


  • (Intercept):1

    • The log-odds of being in the “over a year ago” category relative to “never” are -1.824 when all predictors are at their reference levels.
    • The RRR indicates that the odds of being in the “over a year ago” category are 0.161 times those of “never”.
    • The CIs suggest this estimate is statistically significant.
  • symptd3 or more:1

    • The log-odds of being in the “over a year ago” category relative to “never” increase by 1.127 for those with symptoms 3 or more compared to less than 3.
    • The RRR indicates that the odds are 3.088 times higher.
    • The CIs suggest this estimate is statistically significant.
  • pb:1

    • The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.154 for each unit increase in perceived benefit.
    • The RRR indicates that the odds are 0.857 times lower.
    • The CIs suggest this estimate is statistically significant.
  • hist2yes:1

    • The log-odds of being in the “over a year ago” category relative to “never” increase by 1.063 for those with a family history of breast cancer compared to those without.
    • The RRR indicates that the odds are 2.896 times higher.
    • The CIs suggest this estimate is statistically significant.
  • bse:1

    • The log-odds of being in the “over a year ago” category relative to “never” increase by 0.956 for those with breast self-exam training compared to those without.
    • The RRR indicates that the odds are 2.601 times higher.
    • The CIs suggest this estimate is marginally significant.
  • detcd3 or more:1

    • The log-odds of being in the “over a year ago” category relative to “never” increase by 0.114 for those with detection likelihood 3 or more compared to less than 3.
    • The RRR indicates that the odds are 1.121 times higher.
    • The CIs suggest this estimate is not statistically significant.

  • (Intercept):2

    • The log-odds of being in the “within a year” category relative to “never” are -2.624 when all predictors are at their reference levels.
    • The RRR indicates that the odds of being in the “within a year” category are 0.073 times those of “never”.
    • The CIs suggest this estimate is statistically significant.
  • symptd3 or more:2

    • The log-odds of being in the “within a year” category relative to “never” increase by 2.095 for those with symptoms 3 or more compared to less than 3.
    • The RRR indicates that the odds are 8.123 times higher.
    • The CIs suggest this estimate is highly significant.
  • pb:2

    • The log-odds of being in the “within a year” category relative to “never” decrease by 0.249 for each unit increase in perceived benefit.
    • The RRR indicates that the odds are 0.779 times lower.
    • The CIs suggest this estimate is highly significant.
  • hist2yes:2

    • The log-odds of being in the “within a year” category relative to “never” increase by 1.310 for those with a family history of breast cancer compared to those without.
    • The RRR indicates that the odds are 3.706 times higher.
    • The CIs suggest this estimate is statistically significant.
  • bse:2

    • The log-odds of being in the “within a year” category relative to “never” increase by 1.237 for those with breast self-exam training compared to those without.
    • The RRR indicates that the odds are 3.445 times higher.
    • The CIs suggest this estimate is statistically significant.
  • detcd3 or more:2

    • The log-odds of being in the “within a year” category relative to “never” increase by 0.885 for those with detection likelihood 3 or more compared to less than 3.
    • The RRR indicates that the odds are 2.423 times higher.
    • The CIs suggest this estimate is statistically significant.

4.3 Prediction

4.3.1 Predict the Log Odds

Predict the log odds for the given multinomial logistic regression model (fitmlog1).

# Predicts the log odds and display the first 6 observations
head(predict.vgam(fitmlog1, type = 'link'))
  log(mu[,1]/mu[,3]) log(mu[,2]/mu[,3])
1         -1.2504928         -0.9509763
2         -1.2504928         -0.9509763
3         -0.2411621          0.3053816
4         -1.2504928         -0.9509763
5         -1.2504928         -0.9509763
6         -1.2504928         -0.9509763

  • predict.vgam(fitmlog1, type = 'link'): Predicts the log odds (logit values) for each observation using the model fitmlog1. The type = 'link' option specifies that the predictions should be on the logit (log-odds) scale.
  • head(): Displays the first 6 observations.
  • log(mu[,1]/mu[,3]): The predicted log odds of being in the “over a year ago” category vs. “never”.
  • log(mu[,2]/mu[,3]): The predicted log odds of being in the “within a year” category vs. “never”.

4.3.1.1 Manual Calculation of Log Odds

Manually verify these predictions using the coefficients from the model and the values of the predictors for that observation.

  • (Intercept):1: -1.2504928
  • (Intercept):2: -0.9509763
  • hist2yes:1: 1.0093308
  • hist2yes:2: 1.2565379
head(dat.m1[1:3,])
# A tibble: 3 × 7
    obs    me sympt    pb  hist   bse  detc
  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1     1     0     3     7     0     1     2
2     2     0     2    11     0     1     3
3     3     0     3     8     1     1     3
4.3.1.1.1 1st Observation (hist2 = no)

Logit for “over a year ago” ([1]) vs “never” ([3]):

\[ \text{logit}(P(\text{me3} = 1)) = \beta_0^{[1]} + \beta_{\text{hist2yes}}^{[1]} \cdot \text{hist2} \] \[ = -1.2504928 + 1.0093308 \cdot 0 \] \[ = -1.2504928 + 0 \] \[ = -1.2504928 \]


Logit for “within a year” ([2]) vs “never” ([3]):

\[ \text{logit}(P(\text{me3} = 2)) = \beta_0^{[2]} + \beta_{\text{hist2yes}}^{[2]} \cdot \text{hist2} \] \[ = -0.9509763 + 1.2565379 \cdot 0 \] \[ = -0.9509763 + 0 \] \[ = -0.9509763 \]

4.3.1.1.2 3rd Observation (hist2 = yes)

Logit for “over a year ago” ([1]) vs “never” ([3]):

\[ \text{logit}(P(\text{me3} = 1)) = \beta_0^{[1]} + \beta_{\text{hist2yes}}^{[1]} \cdot \text{hist2} \] \[ = -1.2504928 + 1.0093308 \cdot 1 \] \[ = -1.2504928 + 1.0093308 \] \[ = -0.241162 \]


Logit for “within a year” ([2]) vs “never” ([3]):

\[ \text{logit}(P(\text{me3} = 2)) = \beta_0^{[2]} + \beta_{\text{hist2yes}}^{[2]} \cdot \text{hist2} \] \[ = -0.9509763 + 1.2565379 \cdot 1 \] \[ = -0.9509763 + 1.2565379 \] \[ = 0.3053816 \]

4.3.2 Predict the Probability

Using fitmlog1 model, the log odds for the 1st observation:

  • Log odds for “over a year ago”: -1.2504928
  • Log odds for “within a year”: -0.9509763
# Predicts the probabilities for the first 6 observations
head(predict.vgam(fitmlog1, type = 'response'))
  over.a.year.ago within.a.year     never
1       0.1711957     0.2309783 0.5978261
2       0.1711957     0.2309783 0.5978261
3       0.2500000     0.4318182 0.3181818
4       0.1711957     0.2309783 0.5978261
5       0.1711957     0.2309783 0.5978261
6       0.1711957     0.2309783 0.5978261

  • predict.vgam(fitmlog1, type = 'response'): Predicts the probabilities for each outcome category for the observations using the model fitmlog1. The type = 'response' option specifies that the predictions should be on the probability scale.
  • head(): Displays the first 6 observations.
  • over.a.year.ago: Probability of being in the “over a year ago” category.
  • within.a.year: Probability of being in the “within a year” category.
  • never: Probability of being in the “never” category.

The calculations below demonstrate how to use the predicted log odds to derive the probabilities for each outcome category, ensuring that the probabilities sum to 1 for each observation.

4.3.2.1 Probability of Being in the Reference Group (“Never”)

\[ P(\text{me3} = \text{never}) = \frac{1}{1 + \exp(\text{logit}_{\text{over a year ago}}) + \exp(\text{logit}_{\text{within a year}})} \] \[ P(\text{me3} = \text{never}) = \frac{1}{1 + \exp(-1.2504928) + \exp(-0.9509763)} \] \[ = \frac{1}{1 + 0.2865693 + 0.3867897} \] \[ = \frac{1}{1.673359} \] \[ = 0.5978261 \]

4.3.2.2 Probability of Being in the “Over a Year Ago” Group

\[ P(\text{me3} = \text{over a year ago}) = \frac{\exp(\text{logit}_{\text{over a year ago}})}{1 + \exp(\text{logit}_{\text{over a year ago}}) + \exp(\text{logit}_{\text{within a year}})} \] \[ P(\text{me3} = \text{over a year ago}) = \frac{\exp(-1.2504928)}{1 + \exp(-1.2504928) + \exp(-0.9509763)} \] \[ = \frac{0.2865693}{1 + 0.2865693 + 0.3867897} \] \[ = \frac{0.2865693}{1.673359} \] \[ = 0.1711957 \]

4.3.2.3 Probability of Being in the “Within a Year” Group

\[ P(\text{me3} = \text{within a year}) = \frac{\exp(\text{logit}_{\text{within a year}})}{1 + \exp(\text{logit}_{\text{over a year ago}}) + \exp(\text{logit}_{\text{within a year}})} \] \[ P(\text{me3} = \text{within a year}) = \frac{\exp(-0.9509763)}{1 + \exp(-1.2504928) + \exp(-0.9509763)} \] \[ = \frac{0.3867897}{1 + 0.2865693 + 0.3867897} \] \[ = \frac{0.3867897}{1.673359} \] \[ = 0.2309783 \]

5 Practical 3: Alternative Approach

This exercise demonstrates an alternative approach to fitting a multinomial logistic regression model using the multinom function from the nnet package. It includes steps to fit the model, change the reference category, extract p-values, and compute confidence intervals.

In the VGAM::vglm function, the reference group is the category with the largest number of observations. In contrast, the nnet::multinom function uses the category with the smallest number of observations as the reference group.

5.1 Estimation

5.1.1 Fitting the Model

# Fit the multinomial logistic regression model
mlog_nnet1 <- multinom(me3 ~ hist2, data = dat.m3)
# weights:  9 (4 variable)
initial  value 452.628263 
final  value 396.169969 
converged
# Summary of the model
summary(mlog_nnet1)
Call:
multinom(formula = me3 ~ hist2, data = dat.m3)

Coefficients:
              (Intercept)   hist2yes
within.a.year   0.2995173  0.2470274
never           1.2504941 -1.0093349

Std. Errors:
              (Intercept) hist2yes
within.a.year   0.1662460 0.413737
never           0.1428933 0.427500

Residual Deviance: 792.3399 
AIC: 800.3399 
Note

Note that the reference group here is “over.a.year.ago” instead of “never.”

5.1.2 Changing the Reference Category

# Create new dataset and change the reference category
dat.m4 <- dat.m3
dat.m4 <- dat.m4 %>% mutate(me4 = relevel(me3, ref = "never"))

# Recheck the order of categories
levels(dat.m4$me4)
[1] "never"           "over.a.year.ago" "within.a.year"  
# Fit the model with the new reference category
mlog_nnet2 <- multinom(me4~ hist2, data = dat.m4)
# weights:  9 (4 variable)
initial  value 452.628263 
iter  10 value 396.169969
iter  10 value 396.169969
final  value 396.169969 
converged
# Summary of the new model
summary(mlog_nnet2)
Call:
multinom(formula = me4 ~ hist2, data = dat.m4)

Coefficients:
                (Intercept) hist2yes
over.a.year.ago  -1.2504803 1.009384
within.a.year    -0.9509794 1.256531

Std. Errors:
                (Intercept)  hist2yes
over.a.year.ago   0.1428926 0.4275097
within.a.year     0.1277115 0.3746633

Residual Deviance: 792.3399 
AIC: 800.3399 
Note

The reference category is changed to “never,” which affects the interpretation of the coefficients.

5.2 Inferences

5.2.1 Getting P-Values

# Calculate z-values
z.test <- summary(mlog_nnet2)$coefficients / summary(mlog_nnet2)$standard.errors

# Calculate p-values
p.val <- (1 - pnorm(abs(z.test), 0, 1)) * 2

colnames(p.val) <- c('p-val intercept', 'p-val hist')
p.val
                p-val intercept   p-val hist
over.a.year.ago    0.000000e+00 0.0182218315
within.a.year      9.592327e-14 0.0007972127

  • z.test: z-values calculated by dividing coefficients by their standard errors.
  • p.val: Two-tailed p-values calculated from the z-values.

5.2.2 Confidence Intervals

# Calculate confidence intervals
confint(mlog_nnet2, level = 0.95)
, , over.a.year.ago

                 2.5 %     97.5 %
(Intercept) -1.5305447 -0.9704159
hist2yes     0.1714807  1.8472880

, , within.a.year

                 2.5 %     97.5 %
(Intercept) -1.2012893 -0.7006696
hist2yes     0.5222045  1.9908576

  • Confidence intervals: Provide a range within which the true coefficient values are expected to lie with 95% confidence.

6 Acknowledgment

I would like to express my gratitude to Professor Dr. Kamarul Imran Musa, Medical Epidemiologist and Statistician, and Professor in Epidemiology and Biostatistics at the School of Medical Sciences, Universiti Sains Malaysia, for his teaching and guidance on the subject matter.

7 References

8 R Codes

The R codes used in this analysis are available at the following GitHub repository: DrPH-Epidemiology-Revision. This repository includes all scripts and data necessary to replicate the analyses presented in this work.

 

A dedicated work by Dr Muhammad Saufi